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We propose a new method to study the quasi-normal modes of rotating relativistic stars. Oscil- 
lations are treated as perturbations in the frequency domain of the stationary, axisymmetric back- 
ground describing a rotating star. The perturbed quantities are expanded in circular harmonics, 
and the resulting 2D-equations they satisfy are integrated using spectral methods in the (r, ^)-plane. 
The asymptotic conditions at infinity, needed to find the mode frequencies, are implemented by gen- 
eralizing the standing wave boundary condition commonly used in the non rotating case. As a test, 
the method is applied to find the quasi-normal mode frequencies of a slowly rotating star. 
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CN . I. INTRODUCTION 

^ • Non radial oscillations of compact stars can be excited in several astrophysical events. For instance, after a neutron 
\ star (NS) is formed in a gravitational collapse, or in processes that may occur during its subsequent evolution; these 
OO ■ include starquakes, glitches, interactions with a stellar companion, or phase transitions to quark matter or to a kaon 
T-H and/or pion condensate, that may arise in the inner core of a NS if the density exceeds some critical value. All these 
. phenomena induce perturbations which set the star in oscillation and, according to general relativity, gravitational 
Q ' waves (GWs) are one of the channels through which energy is dissipated. 

In addition, due to rotation some modes may grow unstable through the Chandrasekhar-Friedman-Schutz mccha- 
>' , nism (CFS instabihty) % these instabilities may have important effects on the subsequent evolution of the star, and 
bi), they may be associated to a further emission of GWs, the amount of which would depend on when and whether the 
growing modes are saturated by non- linear couplings or dissipative processes. It is therefore important to know at 
I ' which frequencies a star pulsates emitting gravitational waves, and to study under which conditions the corresponding 
^ I modes become unstable. 

' If one assumes that the star does not rotate, the mode frequencies can easily be computed by solving the equations 
(N| , of stellar perturbations, which have been formulated in the Sixties and further developed in later years j^, 
0^ ' These equations have been integrated for a large variety of equations of state proposed to describe matter in a NS 
! @ • These studies show that the identification of the frequency corresponding to the excitation of a stellar mode (for 
^\ • instance the fundamental mode which is likely to be the most energetic) in a detected gravitational signal, would 
allow to infer interesting information on the composition of the inner core of a NS and on the equation of state of 
: matter at supranuclear densities. 

However, all stars rotate, and our present knowledge of the quasi-normal mode (QNM) spectrum of rotating stars is 
far to be complete. The perturbative approach which works so fine in the non rotating case, when gen eralized to include 
rotation shows a high degree of complexity, even if the star is only slowly rotating [^, 
^ . , major difficulty arises because, when using the standard spherical harmonics decomposition of the perturbed tensors, 
' modes with different harmonic indexes couple, giving rise to an infinite set of dynamical, coupled equations. For this 
reason, in all studies based on this approach simplifying assumptions are introduced: either the couplings between 
oscillations with different values of the harmonic index / are neglected, or Cowling's approximation is used (i.e. 
spacetime perturbations are neg lected) @, i Ud, lill, m, E [3, To our knowledge, the only place where the 
oscillations of a slowly rotating star are studied without making use of any of these restrictive assumptions is in [l6j 
where, however, only r-modes have been considered. 

An alternative approach consists in solving the equations describing a rotating and oscillating star in full general 
relativity, in time domain. However, current studies based on this approach also make use of strong simplifying 
assumptions, or restrict to particular cases. For instance, the Cowling app roximation has been used in several papers 
[l3|; in [l8l| only quasi-radial modes {I = 0) have been considered; in [3) only the neutral mode (zero-frequency 
mode in the rotating frame) has been studied; in (2]| only axisymmetric (m = 0) modes have been analysed, using the 
conformal flatness condition. In [S^l , the frequencies of axisymmetric modes (m = 0) with I — to I — 3 have been 
computed for rapidly rotating relativistic polytropes, using the Cowling approximation; a comparison of the results 
for / = obtained in the Cowling approximation in [13], with those found in full GR in 18], shows that Cowling's 
approximation introduces large errors in the determination of the fundamental mode frequency. 

In this article we develop a general method to find the quasi-normal mode frequencies of a rotating star. We 
perturb Einstein's equations about a stationary, axisymmetric background describing a rotating star. The perturbed 
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quantities are expanded in circular harmonics e""'^. As we are looking for quasi-normal modes, we assume a time 
dependence e~''^', with oj complex. Due to the background symmetry, perturbations with different values of lo and 
m are decoupled; thus, for assigned values of {uj,m), the perturbed equations to solve are a 2D-system of linear, 
differential equations in r and 9. In this paper we do not derive explicitly the perturbed equations in the general case 
of rapidly rotating stars, since they will be studied in a subsequent paper. Our goal here is to describe the general 
method. 

Two are the main ingredients on which our method is based: i) the perturbed equations are integrated using 
spectral methods; ii) the boundary conditions at the center of the star and at radial infinity are implemented by 
suitably generalizing the standing wave approach which has been used to find the QNM frequencies of non-rotating 
stars [3, [j, m, [23|. These points will be described in Section |IT1 To test our approach, in Section Hill we find the 
frequency of the fundamental mode of a slowly rotating, constant density star, as a function of the rotation rate. 
Concluding remarks are in Section [TVl 

II. THE GENERAL METHOD 

The metric which we consider as a background is stationary and axially symmetric. It can be cast in the general 
form 

(ds2)(o) = g^^^dx^'dx" = e-^dt^ + e'f'{d(f> - ujdt)^ + e^'^dxl + e^'dxl (1) 

where ly, 'ip, M2, ^J■3 are functions of the coordinates X2t X3. In llj there is still some gauge freedom, which allows to 
write the metric in a simpler form, like 

= e-'^^'''>dt^ + e^^'^'HdcP ~ w{r,e)dtf + e'^^'^'^'^idr^ + r^de^) (2) 

as in [1^ , or like 

as in [2^ . [2^ . In the following we shall not specify explicitly the gauge, but we will require that some properties 
are satisfied; in particular we require: i) that the spacetime is described by the coordinates {t,r,9,(j)), ii) that ^ 
are Killing vectors associated with stationarity and axisymmetry respectively, and iii) that 9, cj) are polar coordinates 
on spheres, i.e. that the surfaces t — const, r — const are diffeomorphic (but not isomorphic) to 2-spheres. As a 
consequence of these assumptions (that are fulfilled by ^ and ([3])) any tensor field defined on one of these surfaces can 
formally be expanded in tensor spherical harmonics, even if the spacetime is not spherically symmetric; this property 
will be useful in Section [lI A 21 

The metric and fluid velocity perturbations can be considered as tensor fields in this background; they are expanded 
in circular harmonics e"""^, and Fourier-transformed in time. Since perturbations belonging to different m and different 
frequency lo do not couple, in what follows m and lu will be considered as fixed, and perturbed quantities will be 
decomposed as follows 

h^^{t,r,9,(f>) 
Su^{t,r,9,(f>) 

Sp{t,r,9,4>) 

5p{t,r,9,cj,) 

The frequency lo is, in general, complex. 

By fixing the gauge, imposing m^m^ = —1 and assigning an equation of state which relates 5p and 5p, the sixteen 
quantities ft.™^^", <5w™'^, Sp™''^, (5p'"'^ reduce to ten {^f™"(r, 6')}i=i^..._io. A possible gauge choice, which we call 
generalized Regge- Wheeler gauge, is described in Appendix [Xj however, other gauges can be considered (see for 
instance [ISj). For simplicity of notation, to hereafter we shall assume that the quantities i/j""^ are scalars with 
respect to rotation (as they are, indeed, in the generalized Regge- Wheeler gauge, see Appendix]^; however, every 
step of the approach we will describe can be applied also to tensorial quantities by a suitable generalization. 

Einstein's equations, linearized about our stationary axisymmetric background, reduce to a system of partial dif- 
ferential equations (PDE) for the functions i?™" in r and 9. To find the QNM frequencies, for assigned values of m 
and uj, we solve these equations by imposing that all metric functions are regular near the center of the star, that the 
Lagrangian perturbation of the pressure vanishes on the stellar surface, and that the solution at infinity behaves as a 



h'J^^^ir, 6')e™'^e-'"* 
(5<""'(r, 6')e™'^e-''^* 



\m(f) —iu)t 



5p"'^{r,9)e 
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pure outgoing wave. The conditions at the center and at the surface of the star can be fulfilled for every value of uj, 
but the outgoing wave condition at infinity is only consistent with a discrete set of (complex) frequencies {wi}; such 
frequencies are the QNM. 

We will now describe how to implement the boundary condition at infinity, and the numerical approach we use to 
solve the perturbed equations. 



A. The boundary condition at radial infinity 

In this section we shall generalize to the rotating case the standing wave approach [l^, [23| used to find the QNM 
frequencies of non rotating stars. To this purpose, it is useful to remind how this method works. 



1. The standing wave approach for spherical stars 

First of all it is worth stressing that by this approach [l^, [23| one can only determine the QNM frequencies of slowly 
damped modes, like the /-, p-, and r-modes; it cannot be applied to highly damped modes, like stellar w-modes or 
black hole's QNM. 

It is known that the equations describing the perturbations of a non rotating, spherical star can be separated 
by expanding the perturbed tensors in tensorial spherical harmonics; outside the star these equations reduce to 
those describing Schwarzschild perturbations, and they can be reduced to the Regge- Wheeler [28|] and the Zerilli f2§| 
equations, for two suitably defined functions which we both indicate as Z''"(r, w). The two wave equations have the 
following form 



0, 



I > 2 



(5) 



where r, is the usual tortoise coordinate and V{r) is a short range potential. The QNM frequencies are the values of 
the complex frequency lo for which the solutions of Equation (O, found by imposing appropriate boundary conditions 
at the surface of the star, behave as pure outgoing waves at radial infinity, i.e. 



A Irn 



(6) 



The standing wave approach consists in the following. Let us assume that Z^™{r,uj) is an analytic function of the 
complex variable w = cr — i/r, and he loq = (Jq — i/tq the frequency of a QNM, with |1/to| ^ (Tq- In general, at radial 
infinity the solution of Equation ([51) is a superposition of ingoing and outgoing waves, i.e. 



(7) 



If cj = LUf), by definition A\^{u!) = 0; since Z'™ is analytic and since |1/to| <C (Tq, we can expand A\™{a) near the real 
CTo as follows 



1 1m 



(loo) ^ A^\ao) - 
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An '(c^o) , 



where " ' " indicates differentiation with respect to a; then, by imposing A\™{luo) = we find 

4™'(^o)=-iro^^r;K). 
Using this relation the function A'™((t), near cr = (Tq (with a and ctq real), can be written as 

i 



4"(^o) + (^ - ^o)A^ '{ao) = -iroA^iao) 



(cr -ao) + 



To 



from which it follows 



{a - (7o f + ^ 



(8) 



(9) 



(10) 



(11) 



where B is a constant which does not depend on a. Thus, to find the frequencies of the QNM it is sufficient to 
integrate the wave equation ([5]) for real values of the frequency cr, and find the values cr^ for which the amplitude of 
the standing wave (jlip has a minimum: these are the QNM frequencies. The corresponding damping times can be 

found through a quadratic fit of |A'™(cr)|^. 
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2.. The standing wave approach for rotating stars 



Let us now consider a rotating star. As discussed above, with a suitable choice of the gauge the relevant perturbed 
quantities reduce to a set of quantities which behave as scalars with respect to rotation: ^)}i=i,....io ("^ 

complex). They must satisfy a set of PDE, obtained by linearizing Einstein's equation, which can be integrated once 
the values of these quantities are assigned at the center of the star, i.e. on a sphere of radius tq <^ R (hereafter R is 
the stellar radius) 



ijr"(^o,0)-i?o7"(^). 



(12) 



The Hll^^(9) are subject to constraints, which arise from the analytical expansion in powers of r of the perturbed 
equations, from the assumption of regularity of the spacetime as r — > 0, and from the requirement that the Lagrangian 
pressure perturbation must vanish on the surface of the star. These constraints reduce the number of independent 
quantities from the ten Hq1'^{9) to a smaller number, say N, i.e. {H^{9)}j=i^...M ■ Being these quantities defined on 
a sphere r = vq, they can formally be decomposed in spherical harmonics 



/ — |m| 



(13) 



where the expansion is truncated dXl = L. Therefore the independent solutions of the perturbed equations correspond 
to the following set of • [L — |m| + 1] constants 



with 



i = 1, . 



and 



(14) 



Given these constants, the perturbed equations for the functions Hl"-'^{r,9) can be integrated for r > vq. 

In the wave zone, far away from the star, the far field limit expansion of the metric describing a rotating star shows 
that the metric reduces to the Schwarzschild solution (see for instance [30i], Chap. 19). This occurs because terms due 
to rotation decrease faster than the "Schwarzschild-like" components. Therefore, as when dealing with Schwarzschild 
perturbations, in this asymptotic region we can define the gauge invariant Zerilli and Regge- Wheeler functions, 
Z^^(r, w) and Zj^(r, w), in terms of the perturbed metric tensor, expanded in tcnsorial spherical harmonics with 
I > 2. This tensor is found by integrating the equations describing the perturbed spacetime outside the rotating star. 
The well known asymptotic behaviour of Z^"^(r, w) and Z'^^{r,u!) is 
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{r,cj) = A' 
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A 



Im 

RW out 



(15) 



A (complex) frequency wo belongs to a quasi-normal mode if, for an assigned value of to, the following condition is 
satisfied for any I: 



i.e. if the set of 2 • [L — \m\ + 1] constants 

{A^™„„(c.),A^™^,„(c.)} with / = |m|,...,L 



(16) 



(17) 



vanishes. 

It should be stressed that this is a big difference with respect to the non rotating case: in that case each mode 
belongs to a single, assigned value of I, and there is degeneracy in to. 
For each assigned value of to, we define the vectors 



H 



H. 



|m| + l 1 



ry|m| + l m 
-"2 



/ aI™''"(lj) \ 

Zer in \ / 



and A" = 



A 



I m I + i m 



A 



I ml 



RW in 
|m| + 1 m 

RW in 



(c.) 



(18) 
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the dimensionality of which is A'' • [L — |to| + 1] and 2 • [L — |m| + 1], respectively. Since the perturbed equations are 
linear, these vectors are related by the matrix equation 

A"(w) = M™(tj)H™; (19) 

the constants Hj"^ do not depend on iv. The coefficients of the complex matrix M™(w) have to be evaluated by 
integrating the perturbed equations. 

Equation (fT6|) . which identifies the QNM eigenfrequencies, can be written as 



M™(wo)H"' = VH™. (20) 

A discrete set of QNM's exists if the matrix M is square, i.e. if iV = 2. Thus, equation (pO| is equivalent to impose 

det(M™(wo)) = 0. (21) 

As we will show in the Appendices, by counting the number of independent equations in the cases of spherical stars 
and of slowly rotating stars we find indeed N = 2. We expect the same to hold also for rapidly rotating stars. 

Let us now restrict the frequency to the real axis. Furthermore, we normalize the constants H"* so that the solutions 
iJ™°'(r, 0) (and consequently Z'-^^{r,a), Z'^{r,a)) are real; this is always possible, because the perturbed Einstein 
equations in the frequency domain have real coefficients, as long as cr is real (if we assume that the fluid is non 
dissipative, so that the equations in time domain are time-symmetric). Thus, in the wave zone we have 

Z%{r,a) = A^"w.„We-'^^-+^kwo„*(Oe'^'-* e M. (22) 

The ingoing wave ampHtudes, ^z^^m ^'"mvim ^e found, as shown in by evaluating Zzer and Zji^r at 
different values of r,, fitting the two functions as a superposition of sin err, and cos err,. It should be noted that 
although H are real, the quantities A^zlrim ^'"mvin are complex (they satisfy the conditions A^zlrout ^ (^''zlrinT^ 
^''Swout = ("^Wm)*)' ttius M is complex. For a real, the vectors H™, A™ are related by 

A™(a) =M™(a)H™. (23) 

By expanding equation (j2ip about ljq = Uf) — i/ro (with |1/to| ^ ctq) as we did for spherical stars, we find that if 
cr ~ (To 



detM"((7) =detM(CTo)[l-iTo(CT-CTo)] ^ |detM'"(cr)| cx W (cr - ctq)^ + — • (24) 

V To 

Thus, the QNM frequencies are found by evaluating the (complex) matrix M for real values of the frequency cr, finding 
the minima of the modulus of its determinant. The standing wave approach has thus been generalized to rotating 
stars. 

B. Spectral methods for stellar oscillations 

We now summarize the procedure to solve the 2D-perturbed equations using spectral methods. They are, indeed, 
very powerful to solve 2D-differential equations, and particul arly useful to implement boundary conditions. For a 
general discussion on spectral methods we refer the reader to |3l| . 

1. Chebyshev polynomials 
We expand all functions of r in Chebyshev polynomials: 

oc 

/(x) = ^ a„T„(x) with r„(cos(u)) = cos(ri u) n = 0, 1, . . . (25) 

n=0 

which satisfy the orthogonality relations 

Tjn{x)Tn{x)—j^^== = ^ (1 + 5m0)5mn (26) 

-1 ^/l — x'^ ^ 
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(n, m = 0, 1, . . . , I = 1, 2, . . . ). The variable x G [—1, 1] is related to r G [a, b] by the following equation 

2r-b-a 



X 



b- 



a 



e[-i,l]. (27) 



• Integrals on Chebyshev polynomials can be evaluated using the Gaussian quadrature method |32l |: truncating 
the polynomial expansion at n — K we get 



where the collocation points are 



= cos (^i^tlZ^) n = 0,l,...,K. (29) 
• The derivative of a function can be expressed as follows 

fix) =Y.iD^nan)Tm (30) 



with 



Dku = 

Dk-ln = Dk+ln+'2kSkn k = 2,...,K 

Don - l[D2n + 2Sin] ■ (31) 



• Given a function V{x), and a function f{x) with Chebyshev expansion 

K 



/(x) = ^a„r„(a;), (32) 



n=0 



the expansion of V{x)f{x) is 



if 



F(a;)/(x) = J2 bnTnix) , (33) 

n=0 

where &„ = VnmO-m with 

2 - J ^ 

2. Associated Legendre polynomials 
We expand all functions of 6 in the basis of the associated Legendre polynomials 

P'™(2/) = (-1)"(1 - yy-/^—P\y), (35) 

with 

y = cose* e [-1,1] (36) 
and m fixed. They are related to scalar spherical harmonics: 



Y'^{9,4,) - / 2^ + ia "^); p'"»(cosg)e""'^ if m > 
^ ^ y 47r (/ + m)! ^ ^ 

F'™(6»,0) = (-l)'(r'-")* ifm<0. (37) 
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Therefore, expanding a function in circular harmonics e""''^ and in associated Legendre polynomials is equivalent to 
expand it in spherical harmonics. p'™'s are eigenf unctions of the Laplacian operator, 

dl + cot Ode ^ P''" = + 1)P'™ , (38) 

sin 6 J 

and, assuming for simplicity of notation to > 0, have the following asymptotic behaviour near the z-axis (see psp ): 

P'™ - (1 - y^)"/^ = (sin 6*)™ if 6* ~ 0, TT . (39) 

This is the asymptotic behaviour of any function of 9 which is regular on the z-axis. Let us consider a function 
/(r, 6, (/)), regular in r = and on the z-axis; let us expand it in circular harmonics 

/(r,0,0) = ^/™(r,cos0)e™^ (40) 

m 

The regularity of /(r, 9, (j)) near the 2-axis imphes that (see for instance f33|) 

hm ^7-''-"-^) ^finite. (41) 
e=o,7r (sin 6*)™ ^ ' 

Therefore, the f"^{r,cos9) can be expanded in the polynomials {P''™'}i=\m\.... with to fixed: 

oo 

r(r,y)= 5] a'™(r)P'™(y). (42) 

l—\m\ 

Thus, associated Legendre polynomials P'™ (with to fixed) are a complete basis for all functions of 9 with the 
asymptotic behaviour ([55]) . 

In order to apply the Gaussian quadrature method to associated Legendre's polynomials, we notice that the poly- 
nomials 

- il-y^)U ' (43) 

form, for each to, a complete basis, with orthogonality relation 

^ - - , 9 C/ -I- toV 

pl^^y^pl _ y2)^dy = 5„ , . (44) 

I 21 + 1 [I — my. 

The P'™(2/)'s are a particular case of Jacobi's polynomials Ji°''^\y) defined by [s^l 

1 

Ji{y)Ji,{y){l-yr{l + yf ^Sw, (45) 

-1 

where a = (3 = m. Therefore, Gaussian integration takes the form 

/(2/fc)i^""(yfc) 



r' fm^y^pl^^y^dy ^ k ^^^^^^^ 



(46) 



where yk and Wk are the collocation points and weights for the Jacobi polynomials with a = /3 = to. In particular, 
the coefficients of the expansion (|42|) are 

2/ + 1 (?-to)! ^ /(j/fc)P""(yfc) 
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3. Differential equations and boundary conditions 

Let us consider a one-dimensional, first order, differential equation 

Z'{x) + V{x)Z{x) ^ (48) 

with X e [—1, 1]. If we expand Z{x) in the basis of Chebyshev polynomials, Tn{x), truncating the expansion a,t n = K, 
i.e. 

K 



Z(x) = ^a„T„(x), (49) 



n=0 



the differential equation (|48p becomes an algebraic equation: 



K 



J2{Dnk + Vnk)ak^0 {n^O,...,K), (50) 



k=Q 



where Dnk and Vnk are defined in ([3T|l. ([34|) . 

The boundary conditions needed to solve Equation (|48p can be implemented using the so-called r-method: we cut 
the last row of (ISOl). i.e. we set 



K 



Y.{Dnk + Vnk)ak^O {n^O,...,K-l), (51) 



fe=0 



and replace the row with the boundary condition; for instance, if we know that Z{xq) = zq, the last row of the matrix 
equation will be replaced with 



K 



'^Tk{xQ)ak ^ zq . (52) 



k=0 

The differential equation (j48p . with the boundary condition Z{xq) = zq, thus reduces to a matrix equation which can 
be solved by LU decomposition [32]. 

This approach can easily be generalized to higher order differential equations (by replacing more rows for the 
associated boundary conditions), to systems of coupled differential equations, and to partial differential equations in 
r,6. In this case, each function is expanded in the basis {r„(a;), P'™(y)}„^; as follows 

K L 

r{x,y) = E E a5:"r„(a;)P'"(y), (53) 

n=Ol=\m\ 

where sums over Chebyshev's and associated Legendre's polynomials are truncated to the orders K and L, respectively; 
the coefficients aj[" can be organized as the components of a vector, defining the collective index 

n) = [l- |m|] {K + l)+n + l = l,2,...,[L-\m\ + l]{K + 1) (54) 

and setting 



air • (55) 



In terms of the expansion (j53p . PDEs in r, 9 (in our case, the equations which describe the perturbations of rotating 
stars) transform into an algebraic equation. As an example of the double expansion (|53p . in Appendix I A 31 we show 
how to solve the Regge- Wheeler equation as a partial differential equation in r, 9. 

III. A TEST OF THE METHOD: OSCILLATIONS OF SLOWLY ROTATING STARS 

As a test of our method, we have solved the equations which describe the perturbations of a slowly rotating stars. 
Following [131, in this case the metric and fluid velocity of the unperturbed star are 

(ds2)(°) = gl^^dxt'dx" = e-^^-'Ut^ + e^^'^Ur'^ + r^{de^ + sin^ 9d(l)^) - 2uj{ry sin^ 9dtd(f) 
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where uj{r) describes the dragging of the inertial frames, and all quantities are expanded at first order in the angular 
velocity of the star, f2. We introduce a rotation parameter e defined by 

e ^ n/^/M/R^. (57) 

In the metric (|56p spherical symmetry is broken only by the term oj; when Einstein's equations are perturbed about 
()56p . only terms which are linear in u are retained, i.e. we keep terms up to order 0(e). As a consequence, perturbations 
with index I are coupled with perturbations with indexes I ± 1 through terms that are of order 0(e), the analytical 
form of which can explicitly be derived. Therefore, when we transform the perturbed equations using the double 
spectral decomposition described in Section fll B( the resulting algebraic equations have a particularly simple form: 
the relevant matrix is "almost-block-diagonal", each block corresponding to one value of /; the off-diagonal blocks 
couple I ^ I ±1. 

These equations can also be obtained in a different way, i.e. by expanding in Chebyshev polynomials the system of 
ordinary differential equations in r derived by Kojima in [3] (hereafter, Kl). This follows from the fact that Kojima's 
equations are derived by expanding the perturbed Einstein's equations in spherical harmonics. 

The general structure of Kojima's equations is the following: 

CP°' [Hi;" , if'" ; ct] = m£ [hIJ" , if'™ ; a] + [^w" ; (j] 

r'^[Z'S,lv;a] = mAA[Zj,"^;a]-fpW[iif i™,if'±i™;a]. (58) 

Here are operators of order 0{e^), which describe the perturbations of the star in the non rotating case; 

£, Af, 'D^^'> are 0(ej[_) operators, which provide the corrections due to rotation. These equations have been 

integrated numerically in [Sl (see also [!l]) using a very strong simplification: the couplings I <—> I ± 1 were neglected 
(i.e. JF*^^^ and V'^^'^ were set to zero). 

Moreover, the equations were solved iteratively, finding the solution for e = first, and then replacing it in the 
terms £[Hq"^, if'™; cr] and N'[Z'^^; a] of eqs. ((58)) . In this way, the right-hand sides of eqs. ((58| become source terms. 

We stress that with our approach we do not need these simplifications anymore; in particular, we do not need to 
neglect the couplings, because we can handle the mixing among perturbations with different Vs using the spectral 
methods and the generalized standing wave approach. 

It should also be mentioned that, when coupling terms are included in the perturbed equations, in order to have 
a good numerical behaviour of the perturbations near the center of the star we need to use a set of variables (in 
particular the O(e^) terms) different from that used in The equations for the new variables are given explicitly in 
Appendix [B] 



A. Comparison with existing results 

As mentioned above, in [§1 (hereafter K2) the equations of stellar perturbations have been integrated for a slowly 
rotating star, neglecting I -f-^ I ±1 couplings, and the QNM frequencies have been found; to reproduce these results 
we have used the same set of variables as in Kl, the same equation of state (EOS) i.e. the polytropic EOS p = Kp^, 
and we have computed the fundamental mode (/-mode) frequency, u/. 

In K2 the real and imaginary parts of cr^ are fitted as functions of the rotation parameter e as follows: 

a] = al,{l + mea'i) + 0{e^) . (59) 

The value of we find, properly normalized, is plotted versus the stellar compactness, M / R, in Figure [1] a). The 
values are in excellent agreement with the results shown in Figure 1 of K2, for n = 1. 

The correction due to rotation, cr^, is plotted versus M /R in Figure[T]b) for different values of e. We find that for 
e < 10"'^ the corresponding curves are indistinguishable, and coincide with the n — \ curve shown in Figure 1 of K2. 
However, for e > 10^'^ different e correspond to different curves, and the fit ([59]) becomes inaccurate: cr^ is no longer 
a constant, and further corrections to ([55)1 are of order O(e^), as expected theoretically. 

In Figure[2]we plot the imaginary part of the /-mode frequency, CTq, and the corresponding rotational correction, ctj, 
as in Figure[TJ cr^ is plotted only for e > 10~^, because for smaller values our numerical approach becomes inaccurate. 
They agree with the values given in Figure 2 of K2, with differences of order O(e^). 
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FIG. 1: The real part of the /-mode frequency of a non rotating, polytropic star with n = 1, is plotted as a function of the 
stellar compactness M/R (a); the frequency shift due to rotation, a'^ is plotted versus compactness for different values of e (b). 
As in K2, couplings among different Vs axe neglected. 
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n=1 




3.2 - b) 
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FIG. 2: The imaginary part of the /-mode frequency (a) and the corrections due to rotation (b) are plotted as in Figure [T] 

B. Including the couplings 



We now apply our method to solve eqs. (|58p in full, i.e. including couplings among different Vs. The numerical 
implementation of the equations presents a problem: there are terms in the equations which depend on and 

on (^f^) : which strongly diverge on the stellar surface. This divergence is particularly problematic when spectral 

methods are used, but it can be cured through a regularization procedure Such regularization goes beyond the 
scope of the present paper, where we only want to discuss a simple implementation of our approach. Therefore, we 
solve the perturbed equations in the case of a constant density, slowly rotating star, such that the divergent terms 
vanish. 

We consider two background models: A, with M/R — 0.2 and B, with M/R = 0.1. Mass and radius fora^ssigned 
values of the central density are given in Table HI The explicit form of the equations and the boundary conditions in 





P (g/cm^) 


M/Mq 


i?(Km) 


M/R 


A 




1.11 


8.08 


0.2 


B 


10^5 


0.40 


5.75 


0.1 



TABLE I: Parameters of the constant density stellar models A and B we use as a background. 

b) 
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r = and r = R are discussed in Appendix[Bl As mentioned above, when couplings are included the equations derived 
in Kl are very unstable when integrated near the center. For this reason we introduce a new set of variables, which 
satisfy a new set of equations shown in Appendix [Bj the appropriate boundary conditions in r = and r — R are 
also shown. Once we assign the value of the harmonic index m, these equations couple polar and axial perturbations 
with |r7i| < I < L. 

The equations have been integrated for m — ±2. We do not set |r7i| < 2 because in that case dipolar {I — 1) 
perturbations have to be taken into account, which are described by equations different from those we consider in 
this paper. We would like to stress that rotational corrections to mode frequencies with m 7^ are much larger than 
those to mode frequencies with m — 0. In K2, Kojima suggested that, at lowest order in e, the frequency shift is 
proportional to m. This is consistent with the results of our numerical integration: the relative frequency shifts found 
in [2l| , where m = perturbations were studied in full general relativity, are an order of magnitude smaller than the 
relative shifts we find for m ~ ±2. 

If the star does not rotate, for any assigned value of I there is a corresponding /-mode frequency, which is the same 
for all m's. If the star rotates, due to the couplings the /-mode belonging to an assigned I acquires contributions from 
different Vs, and its frequency and damping time change. Furthermore, the degeneracy in m is broken by rotation. 

The real part of the /-mode frequency, i/f = a^/{27r), is plotted as a function of the rotation parameter in Figure [3] 
for the two considered stellar models. The solid line represents the frequency of the I — 2 /-mode of the non-rotating 
star. The dashed lines are the frequencies of the lowest lying fundamental mode of the rotating star, with m — 2 and 
TO — —2, assuming L = 4. Our calculations refer to e < 0.05, since for higher values the slow rotation approximation 
becomes inaccurate and the results cannot be trusted anymore. 
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FIG. 3: The real part of the /-mode frequency, Vf = o-f/(2n), is plotted versus the rotation parameter e. The data refer to 
two models of constant density star (see text). We include couplings up to i = 4. The dashed line refers to the non rotating 
star; the dotted lines refer to the modes m = ±2 of the slowly rotating star. 



Figure[l]shows, in a smaller range of the rotation parameter, the frequency Vf computed by truncating the expansion 
in Z at L = 2 (i.e. without couplings), L = 3 and L — A. It is evident that for slowly rotating stars the contribution 
of the couplings is a small correction, and that there is convergence as L grows. 

The relative frequency shift due to the couplings is well approximated by a quadratic behaviour in e: 

as shown in FigurelHl Therefore, the contribution of the couplings is of order O(e^), as argued by Kojima in K2, and 
a J is well described by a quadratic fit of the form 

af = a^il + mea'j, + e'a'^), (61) 

where we should remind that terms of order O(e^) are of the same order of the terms which we are neglecting in the 
perturbed equations ab initio. This fit is accurate up to e 10~^; for larger values a' and a" are no longer constant, 
as shown in Figure [5] as e grows, a' changes linearly, and the change is negative if to > 0, positive if to < 0, yielding 
in both cases a shift to lower values of the total frequency a. This explains the small asymmetry between i>f{m = 2) 
and i^/(to = —2) shown in Figure [31 Notice that deviations from the fit (pTjl are always of order O(e^), consistently 
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FIG. 4: Details of left panel of Figure O we show the contribution of the different couplings to the /-mode frequency, for model 
A. The different curves are obtained by including in the perturbed equations couplings up to Z = L. 
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FIG. 5: Frequency shift as given by equation (|60|) for model A 



with our approximation scheme. For 0.01 < e < 0.05, the coefficients a" are ~ 0.15 (if to = 2) and ^ 0.20 (if m = —2); 
for e < 0.01 they are too small to be correctly extrapolated with our codes. Finally, the damping time of the /-mode 
is shown in Figure [7] as a function of the rotation rate. 

It is worth stressing that, as the stellar rotation increases, the frequency of the counterrotating (i.e. to — —2) mode 
decreases faster than expected by the simple linear fit ([55]) . Furthermore, Figure [7] shows that the damping time of 
the counterrotating mode increases sharply, even for small rotation rates. This indicates that the CFS instability may 
occur for values of the rotation rate lower than expected by simple linear estimates. 

It should be mentioned that the equations derived for the perturbations of a slowly rotating star in Q and in 
Appendix [B1 are not appropriate to study the r-modes, because the frequency a, which is a dimensionful scale of the 
problem, is of the same order as the "small" parameter e. The shift of the r-mode frequency due to slow rotation in 
a relativistic star has been studied in [l6j . taking into account the couplings between perturbations with different Vs. 



IV. CONCLUDING REMARKS 



In this paper we propose a new approach to find the quasi-normal mode frequencies of rotating relativistic stars. 
We describe the main features of the general method, and test it in the particular case of slowly rotating stars. The 
application of our method to rapidly rotating stars will be discussed in a forthcoming paper. 

We give explicit formulae (whose numerical implementation is straightforward) which allow to transform functions 
of r, 9 in vectors, and systems of coupled differential equations (involving derivatives in r, 9) in algebraic, matrix 
equations. Furthermore, we show that, once the 2D-equations describing stellar perturbations have been solved, the 
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FIG. 6: The coefRcients a', a", plotted as functions of the rotation rate for model A. 
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FIG. 7: Damping time of the /-mode as a function of the rotation rate, for model A. 



frequencies and damping times of slowly damped, quasi-normal modes can be found by looking for the minima of the 
determinant of a properly defined matrix, evaluated as function of real frequency, thus generalizing the standing wave 
approach [H, [11] to rotating stars. 

We have tested the method in the case of slow rotation; the system of partial differential equations from which we 
start is formally the same as in Q (apart from a redefinition of some variables described in Appendix |B]). However, in 
our approach we transform that differential system in a system of algebraic, coupled equations; thus, the advantage 
of our method is that it is much easier to handle the couplings among different values of I, which in correspond 
to couplings between different partial differential equations. For this reason we are able to study the shift of the 
fundamental mode due to rotation, taking into account the / <-^- / ± 1 couplings, to our knowledge for the first time in 
the literature. 

In this paper we show that, as the rotation parameter e increases, the frequency of counterrotating modes decreases 
at a rate higher than linear in e. Furthermore, the corresponding damping time sharply increases, even for small 
rotation rates. This suggests that the CFS instability for a generic mode should occur for values of the rotation rate 
lower than expected by simple linear estimates. This result complements what found in p^ . where the equations of 
stellar perturbations where integrated in full general relativity looking for neutral modes, and it was shown that the 
CFS instability sets in for smaller rotation rates than in Newtonian gravity. 

It should be mentioned that an alternative approach to find the QNM frequencies, based on a characteristic 
formulation of the perturbed equations and a complexification of the radial coordinate, which could be generalized to 
rotating stars, has recently been proposed (3^ . 
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APPENDIX A: GAUGE CHOICE 

In this Appendix we denote with A, B, . . . the coordinates t, r, and with a,b, . . . the coordinates 0, (f>; we define 
7ah = diag(l,sin^ 6). 

1. The generalized Regge- Wheeler gauge 

There are many possible gauge choices to study the perturbations of a stationary, axisymmetric spacetime (see 
for instance [TgjV Here we describe a particular gauge which has the property to reduce, when the background 
becomes spherically symmetric, to the well-known Regge- Wheeler gauge [28[. As discussed in Section we assume a 
background metric in the coordinates {t,r,9,(f>) like ([2]) or ([3]), with ^ Killing vectors. The metric perturbations 
have the form 



We shall show that it is possible to fix the gauge in such a way that the metric perturbation takes the form 



(Al) 



V 





) i?{""(r,6') 









and depends on six quantities 



(A2) 



(A3) 

These quantities behave as scalars with respect to rotations. In order to fix the gauge (|A2[) we impose the following 
conditions 



de sin 6* 



rr(^,o)/i:j,-(r,0) 



Y'"\e,Q)K^^{r, 9) 



VZ. 



These conditions correspond to setting to zero four functions of (r, 9) : 



(A4) 
(A5) 



hZ^{r,9) 

OG 

E 



1 



sin^ei 
d6'sin6' 



hZ^{r,9)=0, 



l—\m\ ' 
oo 



E 



d9 sin t 



Y^rid,0)hir{r,9) 
Y'^{9,0)hlTir,9) 



sm 



-Y^"^{9,0)hr^'^{r,9) 



sm 



-Y'"^{9,0)h^^-ir,9) 



= 0, 



= 0. 



(A6) 



and can be imposed through a diffeomorphism generated by the vector field 

e;.(t,r,0,0)=C™"(r,0)e'^ 

which depends on four functions of (r, 9). 



(A7) 
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The relation between (|A2p and (|A4p is trivial, but to show that (IA5p implies (|A2p is less obvious. If we integrate 
by parts (|A5|) . we find 



d6iy'"(6l,0) 



{s\neh%-{r,e)).e + —-nhA4:{r.e) 



0. 



As it holds for all Z's, the term in square brackets must vanish identically, i.e. 



ITO 



(sin0ft:^-(r, 6)) , + —h^:i:{r, 6) ^ 
' sm f 



(A8) 



(A9) 



therefore we can express the four quantities (ji^^{r,9),h^^{r,6)j in terms of two scalar functions h^^{r,6) such 
that: 



'ir,9) 



(AlO) 



This, together with (|A4p . gives (|A2[) . This gauge choice is implicit in the formulation used in Kl to describe the 
perturbations of slowly rotating stars. It can, in principle, also be chosen to describe perturbations of rapidly rotating 
stars. 

We stress that the existence of the generalized Regge- Wheeler gauge, in which all perturbations are expressed in 
terms of functions that are scalar with respect to rotation, is important because it provides a solid basis to the approach 
described in this paper. Indeed, it guarantees that expanding the perturbations in tensorial spherical harmonics is 
equivalent to expand in circular harmonics e'"*"^ firstly, and then to expand the functions appearing in the resulting 
equations in {r,9), in associate Legendre polynomials, P'™(6'). 



2. Relations with the Regge- Wheeler gauge 

In order to better understand how the gauge (|A2p is related to the Regge- Wheeler (RW) gauge, we now expand 
hfivit, r, 9, 0) in tensor spherical harmonics. This is always possible (on a surface t = const., r = const.), but typically 
it is not useful if the background is non-spherical, since the dynamical equations couple perturbations with different 
Vs. Anyway, for slowly rotating stars the couplings are small, and the spherical harmonics expansion, as described in 
Appendix [B] and in [3] , turns out to be useful. 

By expanding the perturbed metric tensor in tensor spherical harmonics, before any gauge fixing we find 





f e''iJ^"(r)y'™(6', 


i>) i7{™(r)r''"(6l,(/)) 
e^iJ^"(r)y'"(6',(/>) 




) + h'^{r)Si"'{e,4>) 
) + /if(r)^i™(0,0) 






K^"'{r)r^'yabY^"'{d, (j)) -f G'™ 


(r)Zi'^(0,0) + /.i™(r)5i'^(0,0) 



(All) 

where 

s'ne, 0) = (5^, 5^™) = (^--i^y;,™, sin^y'r) (A12) 

are the axial vector harmonics, and Zab and Sab are tensor harmonics satisfying ^""^ Zab — "l°'^Sab — 0, with polar and 
axial parity, respectively. 

The RW-gauge for a spherical background, imposes [281] : 

Co;=/»i'"po; = C-G'™=0. (A13) 

If we consider (|A13p in the case of a non-spherical background, we see that it is equivalent to the gauge (|A2p . Indeed, 
expanding liAait, r, 6, (j)) in vector spherical harmonics, we find 

h7a{r. Oy^^ = E [hXoi{r)Y^r{0. </>) + ^a" <^)] , (AM) 

I 
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which, if h}^^^i = 0, reduces to 



mi 



/iA,sin6'/iAe e'™*", (A15) 
smf ' J 

where we have defined the scalar functions 

3. The Regge- Wheeler equation for a spherical star as a 2D-equation in r and 9 

We conclude this section with a simple exercise. Choosing the gauge (jA2l) , we derive the equations which describe 
the axial perturbations of a non rotating star, i.e. the equations for /i™"(r, 6*) (in this case fluid perturbations are 
decoupled from metric perturbations), and we show how to solve them using spectral methods. We shall follow the 
lines of the well-known derivation of the Regge- Wheeler equation [28], with one difference: we do not expand the 
perturbations in spherical harmonics. By defining the function 

Z""^{r, 0) = 6')e(''-^)/2 , (A17) 

from Einstein's equations we find 

/iS""(r, e) ^ ^irZ"'^{r, 0))'e("-^)/2 ^ (^18) 

1(7 

where "' indicates differentiation with respect to r. Z™'^{r,9) satisfies the partial differential equation 

d^^Z""^ + {a^ -V)Z''"^ ^0 (A19) 

where the coordinate r* is defined by 

^ = e(^-'')/2 . (A20) 
dr 



In this formulation, is a differential operator: 

(A21) 



dl + cot edg ^ + Ati{p - p)r^ 

sin 9 J r 



The usual one-dimensional Regge- Wheeler equation can easily be recovered by expanding Z™ (r, 9) in scalar spherical 
harmonics: 

^'^ \v^ = ^-{l{l + l)-'-^+An{p-py) ^ ^ 

where 

Imin = max(|m|, 2) , (A23) 

and Z^™'^{r) is the standard Regge- Wheeler function. 

Let us briefly describe how equation (|A19p can be integrated using spectral methods and the standing wave approach. 
We shall integrate this equation for real values of the frequency, therefore in the following we shall set lo ^ a. The 
integration range in r* is = [r*i,r*2], i.e. given Z™'^(r*i,6') we want to know Z™''{r^,2,9). The starting point r*i 
corresponds to a small sphere near the center of the star, with r = ro <C i?, where the Regge- Wheeler function is given 
by an analytical expansion (see below). The final point, r^2, corresponds to a point in the wave zone, r = Too S> R, 
where the ingoing and outgoing amplitudes can be extracted. 
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We rescale the variables r» and as follows: 

2r, - r,i - r»2 



cos6ie [-1,1], (A24) 



so that r*i corresponds to a; = 1, and r*2 to a: = — 1. Then, we perform the double expansion ([53|) of the Regge- Wheeler 
function Z™'^(r, 0) defined in (|Al7p . for assigned values of m, a: 



K 



Z-^%x,y) = Y, a':^-T„ix)P'^{y). (A25) 



n=0 l=l„ 



The expansions in Chebyshev's and associated Legendre's polynomials are truncated at K and L, respectively (for 
instance, L = 10, K — 20). The boundary conditions at the center of the star are imposed by assigning ^'"'^ and its 
derivative at a; = 1: 



L 

dr,Z""'{x^l,y) = Y ^H^'^P^^'iy), (A26) 



where the analytic expansion of the Regge- Wheeler equation gives [3] 



'0 — ^0 ^^m\'^ro 



1+21 



z[ = e-(0) [{I + l)r[, + + 3)XL(rT)ro ^ 

X' r.^ - a + 2) [i(2?-iM0)-p(0)]-aV(") 

^-^'"^ = 2(27T3^ ^^^^^ 



(note that while in Zq ^ and Xj-„, the index I is a superscript, in Tq it is an exponent). 
The constants -ff'™ form a vector 



jjm = / ^(m ^ (A28) 



which can be freely assigned; for each vector H we have one solution of the equation. Notice that we have one constant 
jijim value of / i.e., in the language of Section Hi A 2[ iV = 1 for the axial parity perturbations. If, in addition 

to axial perturbations, polar parity perturbations are considered, described outside the star by the Zerilli equation, 
then there is another constant to be assigned for each value of I. Therefore, if all metric perturbations are considered, 
TV = 2 as discussed in Section Hi A 21 

We now project equation (|A19p in the basis of Chebyshev's and associate Legendre's polynomials. The operator 



dl + cot 9de 5- (A29) 

sin 9 



is diagonal in the P'*" basis, with eigenvalue + Therefore, in this basis the operator V defined in (|A21[) reduces 
to the one-dimensional Regge- Wheeler potential V^{r) (|A22|) . We introduce (as in in Section fllB 11 Equations (|3ip . 
([34]) ) the derivative matrix D„„' and the potential matrix Vnm obtained projecting y'(r) on Chcbyshev polynomials. 
To write the matrix equation we define the collective index 

i(l,n)^{l-\m\){K + l)+n + le (A30) 

with M = (i — \m\ + 1){K + 1), and we define the A/'-dimensional vector of components 



(A31) 



The matrix equation has the block-diagonal form 

Lu-a'^ = (A32) 
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with 



L 







l—\m\ 



Tin' ' ^ ^ 



i=|m| + l 



(A33) 



l=L 



The boundary conditions are implemented by replacing the last two lines of each block, i.e. each ?(/, K)-i\\ and 
K — l)-th lines, with the conditions at the center in terms of the vector H™ — (i?'™): 



K 



/ Trim 



O.K 



(A34) 



n,k 



By inverting the matrix equation (jA32p we find the coefficients a™ and then 1, y). 

In this way, for each choice of (to, cr), and of the vector of initial conditions, we can integrate from tq to 
The ingoing amplitude at infinity 



A"" {(7,0) ^yl'"'((7)P'"'(6') 



(A35) 



can be found using the algorithm described in |24] . Therefore, for each choice of the vector of initial conditions 
H™ — (iJ'™), we have the vector of ingoing amplitudes at infinity 



H™ = (0,0,...,1) 



A™ = (y^'™(cr)). 

This procedure is linear at any step, thus, repeating it for all H™: 

H™-(1,0,...,0), H™ = (0,1,...,0), ... 
we find the expression of the matrix M™ = (yVJ™!" (cr)) defined by 

yt'™(o-) Al"l"'(CT)ij'''". 

As explained in Section [II A 21 near a mode a ^ ao the modulus of the determinant of this matrix behaves as 



|detM'"(a)| cx \ (a ~ aoY 



To 



(A36) 



(A37) 



(A38) 



(A39) 



and the frequency (Tq and the damping time tq of the mode can be found by a quadratic fit in a. 



APPENDIX B: EQUATIONS FOR THE PERTURBATIONS OF SLOWLY ROTATING STARS 

This derivation of the equations describing the perturbations of a slowly rotating star is based on the work of 
Kojima [7|, denoted as Kl. As discussed above, since we are considering slowly rotating stars, we can first expand 
the perturbations in spherical harmonics, getting a system of coupled ODE in r, and then expand this system in 
Chebyshev polynomials, getting an algebraic matrix equation. The first step of this program is equivalent to the 
derivation of Kl, with one difference: we are going to reformulate the equations in terms of a different set of variables, 
which are numerically well behaved near the center of the star. Following Kl, we shall assume I > 2. 

The background configuration, describing a slowly rotating, stationary and axially symmetric star is given by 
equation (|56p . The pressure p and the energy density p are found by solving the TOV equations; we assume that the 
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equation of state of matter in the star is barotropic, p — p{p)\ therefore, = = ^- The perturbations of the 

background (j56p can be written as 



1 771 

5u^ = 







,) /i[,'"(r)^^'"(0,0) 






e^H!p{r)Y^'^{e,(j)) 














,0) 
Ki^{r)r^ sin^ 6'r'™(6', 







-A 








-A 


47^(p^ 






-A 



(T/''"(r)Y;^"(6',0) + J7'™(r)S'^™(6',,/.)) e-'""' 



rim ( 



6p = 5p^"'{r)Y^" 



(Bl) 



with a real. The hnearized Einstein equations for the radial part of these quantities are ordinary differential equations 
in r. As explained in Section [lIH the general structure of these equations is: 



(B2) 



The quantities i/;^™, i?'", V'™, J/'", (Jp'", can be expressed in terms of i?^™, X'™, once equations 

(IB2I) have been solved. 



1. The 0(e°) equations 
The equations at 0{e^) describe the perturbations of a non-rotating star. 

The equations for polar perturbations inside the star are a system of two second order ODE in X'™, iJp™: 

pA 

+4^[3Mr - 47rpr* - e^(Af + \'Kpr^f\H^^ = (B3) 



pA -2 _ 1 



+ ^[(nr + 4M + 87rpr3)c;2 - (n + 2)r + STrpr^Ji/^™ = . (B4) 



Once -ffQ™, ii''™ have been determined, and i?2™ can be computed through the following relations: 

iji" = -—(h\"''- if'™' + ^(Af + Anpr^)Hl." 
i(T \ 

iJ^" = iJ^'". (B5) 
Other relations (see Kl) give the fluid perturbations i?'"*, V''"^, i7'™, (5p'™, (5^'™ in terms of the metric perturbations. 
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Outside the star, the polar perturbations reduce to a system of three first order ODE in K^"'^, H[™ and Hjf^: 



- {l}lm 



- Tjlni 



r r 



2Me^ 



2 A " + 1 

a re 



i7 ™ = , 



-m 



where we have defined 



Im 



in order to have equations with real coefficients. Furthermore, there is an algebraic constraint: 



-(4)/mr 



{nr + 3M)rH^"' - [a^r'^ - [n + l)Mr]H'{^ = 



(B6) 
(B7) 

(B8) 

(B9) 



(BIO) 



Equations (|B6P - (jB8P are equivalent to the Zerilli equation, but they have the advantage to be easily generalizable to 
rotating stars, as we will see below. The Zerilli function Zzer can be computed in terms of K, Hi : 



rvlm 



Ira V fjlm 



nr + 3M 

The axial perturbations are described by the Regge- Wheeler equation 



^ [1(1 + I) 6M ^ . 



ryira ri 

^RW — ^ ' 



where the coordinate r* has been defined in (jA20p and 
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{Z^Rwr)' 



(Bll) 



(B12) 



(B13) 



An analytical expansion of equations (|B3p . (|B4p . (|B12p near the center of the star gives, for each value of Z, three 
independent conditions at the center: 
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^1+2 ^ f^^lm(+)^l+4 
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(B14) 



H, 



■'RW 



^1+2 ^ J,l,lm(-)^l+i 





(B15) 



rylra 
^RW 








(B16) 



where the expressions for kh^"^^^\ k^"^'^^\ z'™ can be evaluated from the analytical expansion. 

We notice that if'™ and Hq"^ behave, as r ^ 0, like r', while the combination iiT'™ — Hq™' behaves as r'+^. 
Consequently, when we expand K''"^ and H^"^ in powers of r about r = 0, we find that the leading terms are 
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coincident. In other words, the differential equations for the variables K , Hq and are linearly dependent 

near the origin (see the discussion in 0]). To avoid this problem, we use as integration variables K^"^ — Hq"^, K^™^ 
and Z^^. 

Finally, we notice that at order 0(e°), i.e. for a non rotating star, equation (|B5P establishes that = Hq^. 
Therefore it is equivalent to use either ifg™ or iJj™- 



2. The O(e^) equations 



We define 



1(1+3) 
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m 



{l-rn){l+m) 
(2/-l)(2;+l) 



^ (;-2)(i+i) 



Ql+1 m 



{l+l-m)(l+l+m) 
(2; + l)(2;+3) 



(B17) 



The perturbed equations inside the star have the form 



[Ho,K] = - 



(n+ l)cr 



+e 



-(A+i/)/2 



iQ 



l — l m 



RW 



L ^"'[Zrw] 



N'"'[ZRw]+e 



iQ 



cr(n — n+) 



(A+3i/)/2 



a{n — n^) 



(B18) 



(J = 1, 2), where "[ffo, i^], are the operators defined in equations (|B3]), jll}, ((BT2)) and E'--'^^"'[Ho, K], 

£)(,/)i±imj^^^j^ 7V'™[Zfl^y], _F'*^™[i7o, are operators at first order in e, whose explicit expressions are given in 
Kl. The operators D^'^^^^^'^lZnw], F''^^"^[Ho, K] couple perturbations belonging to different Ts, and were neglected 
in the numerical integration of K2. 

In order to have equations with real coefficients, we need to get rid of the factors i in l|B18[) . To this purpose, we 
redefine the Regge- Wheeler function by a factor — i 



AZ 



Im 

RW ' 



thus equations (jBlSp become 



(B19) 



- (J)/m 



(n + l)cr 



+e 



'(A+i/)/2 
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L™'[Zrw] 
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a(n — n_) 

,(A+3<.)/2 



a{n — n^) 



am — n_ 



(j(n — n+) 



(B20) 



This rescaling consistently eliminates all imaginary units from the equations. 

The numerical integration of equations (jB20p presents a serious problem. If we perform an analytical expansion 

near the center of (lB20l) . we find that the coupling terms D^-^'^^-^"'[Z%,] become larger than L,^;[f "[iJ^™, if'™] as 
r 0. The reason behind this pathological behaviour is that near the center of the star 



K' 



Im 



K 
Hi 



^1+2 + ( (9(g) terms ) • 



Ira 



(B21) 



If we use as integration variable iJj™ instead of H^"' (notice that while in the non rotating case H^"' — Hi^^, if the 

0(e)) this problem is overcome, since 



star rotates Hq"^ = iJj 



K 



Im 



Tjlm 
^2 



J+2 



(B22) 
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Consequently, the coupling terms in the perturbed equations are smaller than L^^f"^[H2",K^'^]. For this reason, we 
have expressed our equations inside the star in terms of K'™ — Hlf^, if'™, Z^^^: 
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''int 
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(n+ 1)(T 



+e 
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-(A+i^)/2 
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(j(n — n_) 
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(B23) 



The operators ^('')'™, £)(-')'±i'» are different from £;("')'™, i:)(J);±i'n 

given in Kl. Their expressions are the following: 



2" 
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+ia{n + l)e^C(3)'™(r) - ia(n + l)C(2)'™(r) 



(B24) 
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-iae(^+'')/2 f- + A' - 2J.') [-2(n - 2n± - 2)x^'^'^' {r) 
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(n - n±)i(2)'±^"(r) + (n - n±)(n - n± - 1)B(2)«±i 



(B26) 
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(B27) 



where /, g, C^"') , a^"') , /3^-'^ , r/("^) , C^'''' , C^'^ , A^^^ , B^^^ are quantities which depend on the perturbations iJl™, K^"", etc., 
and which are given in Appendix B of Kl. 
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At the surface of the star we compute Hq^ and the other perturbations in terms of H!f^, K^"^, Z^^. We impose 
the vanishing of the Lagrangian pressure perturbation (see Q). This reduces the number of freely assigned constants 
from three (times L — \m\ + 1), which correspond to the three independent solutions (jB14P - ()B16|) . to two (times 
L — |m| + 1), i.e. iV = 2 as discussed in Section Hi A 21 

Finally, the equations in vacuum are, as in Kl, 

a \ n — n_ n — / 

L^m^Znw] = - (mN^"^[Znw] - - Ql+lJlli^ll^llE^] (B29) 

a \ n — 71- n~ n+ I 

(J = 1, . . . , 4) where 1^1''"' [Hq, Hi,K], L'"^[Zrw] are the operators defined in (|B3|1 - ([BT0)) . ((BT2)) . and the expressions 
of E^J)h''[Ho,K], M^^^^^IZrw], F^-'^^"^[Ho,K] are given in Kl. 

When r ^ R, the background spacetime is with good approximation spherically symmetric, because the terms due 
to rotation decrease faster than the "Schwarzschild-like" components (see for instance [sOl, Chap. 19). Therefore, 
spacetime perturbations satisfy the Zerilli and the Regge- Wheeler equations. 

In this limit, equation (|B29p becomes the Regge- Wheeler equation for the function whereas the Zerilli function 

is related to the solution of equation (|B28p by (jBlip . At radial infinity, the amplitude of the stationary wave 
(^Zer in (^) ' ^flw in ( g)) C an be computcd in terms of and . We can then apply the stationary wave approach 
described in Section fll Al and in Appendix lA 31 

Equations (|B23|) . (jB28p . (|B29p can be integrated using the spectral decomposition in Chebyshev's polynomials 
as explained in Section IIIBI There is a main difference with respect to the example described in Appendix lA 3[ 
which refers to the axial equation for a non rotating star. While the matrix (jA33P is block-diagonal - each block 
corresponding to a value of I - the matrix representing equations (jB23p . (jB28p . (jB29p presents, in addition to the 
block-diagonal terms, components of order 0(e), which couple ^ <-> Z ± 1. 

All equations in this paper have been checked using Maple, and we have made several cross checks in order to be 
sure that the Fortran implementation of the long expressions (jB24l) - (jB27p are correct. 
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